#version 120
#extension GL_EXT_gpu_shader4 : enable //access to integer and bit operation
uniform mat4 invViewMatrix;

//uniform sampler3D volumetexture;
uniform usampler3D nodePool3D;
uniform sampler3D blockPool;

uniform sampler1D tranf;

uniform float stepsize;
uniform float begin;
uniform float isoValue;
uniform float var;



vec4 DirectVolumeRendering(vec3 front, vec3 back, float upperThreshold, float threshold, float offset)
{
	vec3 dir;
	dir = back - front;
	float len = length(dir);
	vec3 norm_dir = normalize(dir);
	//float iso_value = 0.1;

	vec3 rayIncrement = stepsize*norm_dir;

	vec3 pos = front;
	vec4 value=vec4(0,0,0,0);
	// accumlate value
	vec4 result = vec4(0,0,0,0);//RGBA
	float raylength = 0;

    float res = 0;
    
	while(result.a < 1.0 && raylength < len)
	{
		float factor = 0.5;

	
		value = texture3D(blockPool, front+raylength*norm_dir);
		raylength += stepsize;

		vec4 texValue = texture1D(tranf, value.a);

		result.rgb += (1.0-result.a)*texValue.a*texValue.rgb;
		result.a += (1.0-result.a)*texValue.a;

	}
	return result;
	//return vec4(front,1);

}

vec2 uv = (gl_TexCoord[0].xy/gl_TexCoord[0].w);
//=====Parameter
const float voxelMax = 255.0;
//1G
const uint nodePoolDim = 28u;
//const uint blockPoolDim = 23u;
const uint blockPoolDim = 25u;
const float invBlockPoolDim = 1.0/float(blockPoolDim);
const uint blockSize = 33u;
//============
struct nodeTexel
{
	uint first;
	uint second;
};


//Utility
nodeTexel intTexture3D(usampler3D tex, uint z, uint y, uint x, uint maxZ, uint maxY, uint maxX)
{
	//maxZ (1:maxZ) number of element = maxZ
	// z = (0: max-1)
	float indexX = (float(x)+0.5/float(maxX))/float(maxX);
	float indexY = (float(y)+0.5/float(maxY))/float(maxY);
	float indexZ = (float(z)+0.5/float(maxZ))/float(maxZ);//TODO: convert to a vec operation

	uvec4 uSample = texture3D(tex, vec3(indexX, indexY, indexZ));

	nodeTexel node;
	node.first = uSample.r;
	node.second = uSample.a;
	return node;
}


uvec3 getNodeClusterPointer( nodeTexel node)
{
	uvec3 cPointer=uvec3(0u,0u,0u);
	cPointer[0] = node.first & uint(0x3FF) ;
	cPointer[1] = (node.first>>10u) & uint(0x3FF);
	cPointer[2] = (node.first>>20u) & uint(0x3FF);
	return cPointer;
}

nodeTexel getChildrenNode(uvec3 clusterPointer, uint k, uint j, uint i)
{
	nodeTexel node = intTexture3D(nodePool3D, clusterPointer[2]*2u+k, clusterPointer[1]*2u+j, clusterPointer[0]*2u+i,
													nodePoolDim,nodePoolDim,nodePoolDim);

	return node;
}

uvec3 getBlockPointer( nodeTexel node) 
{
	uvec3 bPointer=uvec3(0u,0u,0u);
	bPointer[0] = node.second & uint(0x3FF) ;
	bPointer[1] = uint( (node.second >> 10u) & uint(0x3FF));
	bPointer[2] = uint( (node.second >> 20u) & uint(0x3FF));
	return bPointer;
}

bool hasBlock( nodeTexel node)
{
	bool block = ((node.first>>30u &  uint(0x1)) >  uint(0))? true:false;
	return block;
}

bool isSubDividable( nodeTexel node)
{
	bool dividable =  ((node.first>>31u &  uint(0x1)) >  uint(0) )? true:false;
	return dividable;
}


struct ray
{
	vec3 o;
	vec3 d;
};

struct range
{
	float tNear;
	float tFar;
};

struct BlockBox
{
	vec3 boxMin;
	vec3 boxMax;
	bool hasBlock;
	uvec3 blockPointer;
	float sampleRate;
	float singleColor;
};

range intersectBox(vec3 entryPoint , vec3 rayDir, vec3 boxmin, vec3 boxmax)//, float *tnear, float *tfar)
{
    // compute intersection of ray with all six bbox planes
    vec3 invR = vec3(1.0) / rayDir;
    vec3 tbot = invR * (boxmin - entryPoint);
    vec3 ttop = invR * (boxmax - entryPoint);

    // re-order intersections to find smallest and largest on each axis
    vec3 tmin = min(ttop, tbot);
    vec3 tmax = max(ttop, tbot);

    // find the largest tmin and the smallest tmax
    float largest_tmin = max(max(tmin.x, tmin.y), tmin.z);
    float smallest_tmax = min(min(tmax.x, tmax.y), tmax.z);

	range tMinMax;
	tMinMax.tNear = largest_tmin;
	tMinMax.tFar = smallest_tmax;
	return tMinMax;
}

BlockBox findCurrentBlockAndBoundingBox(vec3 pos)
{
	vec3 boxmin = vec3(0,0,0);
	vec3 boxmax;// = vec3(1,1,1);
	float boxsize = 1.0;//vec3(1,1,1);

	uvec3 nodeClusterPointer=uvec3(0u,0u,0u);
	uvec3 childrenIndex;

	nodeTexel currentNode;
	currentNode.first = 1u<<31u; currentNode.second = 0u;

	
	float scaleFactor = 1.0;
	
	bool SubD=true;

	BlockBox renderBlock;
	renderBlock.sampleRate = stepsize;

	while(SubD) //@@can't write isSubDividable(currentNode);@@
	{
		//decending
		childrenIndex.x = uint((2*pos.x));
		childrenIndex.y = uint((2*pos.y));
		childrenIndex.z = uint((2*pos.z));

		boxsize *= 0.5;
		boxmin += vec3(childrenIndex)*boxsize;
		//scaleFactor *= 0.5;

		pos = fract(2*pos);
		
		currentNode=getChildrenNode( getNodeClusterPointer(currentNode), childrenIndex.z, childrenIndex.y, childrenIndex.x);
		SubD = isSubDividable( currentNode );

		renderBlock.sampleRate *= 0.5;

		//scaleFactor *= 0.5;
	}
	uvec3 blockPointer=getBlockPointer(currentNode);
	boxmax = boxmin+vec3(boxsize);

	renderBlock.singleColor = float(currentNode.second)/voxelMax;
	renderBlock.hasBlock=hasBlock(currentNode);
	renderBlock.boxMin = boxmin;
	renderBlock.boxMax = boxmax;
	renderBlock.blockPointer = blockPointer;
	
	return renderBlock;
}

BlockBox findCurrentBlockAndBoundingBox_MultiRes(vec3 pos, float vt, float t, vec3 O) 
//vt is the distance between eye to the image plane, and t is the distance between eye to the current Position
{
	vec3 boxmin = vec3(0,0,0);
	vec3 boxmax;// = vec3(1,1,1);
	float boxsize = 1.0;//vec3(1,1,1);

	uvec3 nodeClusterPointer=uvec3(0u,0u,0u);
	uvec3 childrenIndex;

	nodeTexel currentNode;
	currentNode.first = 1u<<31u; currentNode.second = 0u;

	
	float scaleFactor = 1.0;
	
	bool SubD = true;
	bool isResOK = false;
	float voxelSize = 1.0/float(blockSize);
	BlockBox renderBlock;
	renderBlock.sampleRate = stepsize*32;

	while(SubD && !isResOK) //@@can't write isSubDividable(currentNode);@@
	{
		//decending
		childrenIndex.x = uint((2*pos.x));
		childrenIndex.y = uint((2*pos.y));
		childrenIndex.z = uint((2*pos.z));

		//check if current resolution is enough (the projected voxel size is smaller than a pixel in image space)
		if(1.0/length(boxmin - O)*voxelSize</*thresholdScreenPixelSize*/var)
			isResOK = true;
		voxelSize *= 0.5;

		boxsize *= 0.5;
		boxmin += vec3(childrenIndex)*boxsize;

		pos = fract(2*pos);
		
		currentNode=getChildrenNode( getNodeClusterPointer(currentNode), childrenIndex.z, childrenIndex.y, childrenIndex.x);
		SubD = isSubDividable( currentNode );

		renderBlock.sampleRate *= 0.5;
	}

	uvec3 blockPointer=getBlockPointer(currentNode);
	boxmax = boxmin+vec3(boxsize);

	renderBlock.singleColor = float(currentNode.second)/voxelMax;
	renderBlock.hasBlock=hasBlock(currentNode);
	renderBlock.boxMin = boxmin;
	renderBlock.boxMax = boxmax;
	renderBlock.blockPointer = blockPointer;
	
	return renderBlock;
}
	////calculate blockoffset in the pool
	//uvec3 blockPointer=getBlockPointer(currentNode);

	//vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	//if(!hasBlock(currentNode))
	//	return vec3(0.0,0.0,0.0);

	////boundary shrink factor
	//float shrinkFactor = (float(blockSize)-1)/(float(blockSize));
	//vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;
	////scale the pos
	//return offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset);

//vec4 accumulate(vec3 rayO, vec3 rayDir, vec4 color, float tStart, float tEnd, float sampleRate)

//ray and dir is alway object coodinate
vec4 accumulateSingleColorBlock(vec4 color,BlockBox renderBlock, range r)
{
	//float invboxsize =1.0/( renderBlock.boxMax.r - renderBlock.boxMin.r);
	//float shrinkFactor = (float(blockSize)-1)/(float(blockSize))*invboxsize*invBlockPoolDim;
	////@@Need 
	//vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;
	float tStart = r.tNear;
	float tEnd =  tStart+(r.tFar-r.tNear);
	//accumlate
	float  t=tStart;
	vec4 value, texValue;
	vec3 pos;	
	texValue = texture1D(tranf, renderBlock.singleColor);

	//TODO analytical integrate
	//safegard
	int i=0;
	while(color.a < 1.0 && t < tEnd)
	{
		if(i>200)
			break;
		i++;
		//pos = blockO+t*rayDir;
		//value = texture3D(blockPool, offset + pos*shrinkFactor + shrinkOffset);
		
		color.rgb += (1.0-color.a)*texValue.a*texValue.rgb;
		color.a += (1.0-color.a)*texValue.a;
		t += stepsize;

	}

	color.rgb += (1.0-color.a)*texValue.a*texValue.rgb;
	color.a += (1.0-color.a)*texValue.a;

	return color;

}


vec4 accumulateBlock(vec3 rayO, vec3 rayDir, vec4 color,BlockBox renderBlock, range r)
{
	//warping the coordinate
	uvec3 blockPointer = renderBlock.blockPointer;

	float invboxsize =1.0/( renderBlock.boxMax.r - renderBlock.boxMin.r);
	//vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	//float shrinkFactor = (float(blockSize)-1)/(float(blockSize));
	//vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;

	vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	float shrinkFactor = (float(blockSize)-1)/(float(blockSize));
	vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;


	if(r.tNear>r.tFar)
		return vec4(0,0,0,1);

	//accumlate
	float tEnd =  r.tFar;
	float  t=r.tNear;
	//float t+=residual;
	vec4 value, texValue;
	vec3 pos;	

	int i=0;
	while(color.a < 0.99 && t < tEnd)
	{
		//safeguard
		i++;
		if(i>200)
			break;

		pos = (rayO+t*rayDir);
		//return vec4(pos,1);
		
		pos = (pos-renderBlock.boxMin)/0.03125;
		//for(int f=0; f<5; f++)
		//	pos = fract(2*pos);

		//return offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset);
		
		value = texture3D(blockPool, offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset));
		//value = texture3D(blockPool, offset + pos*shrinkFactor*invBlockPoolDim + shrinkOffset);
		texValue = texture1D(tranf, value.a);

		color.rgb += (1.0-color.a)*texValue.a*texValue.rgb;
		color.a += (1.0-color.a)*texValue.a;
		t += stepsize;

	}
	//color.rgb += (1.0-color.a)*texValue.a*texValue.rgb*(tEnd-(t - stepsize));

	//visualize block Pointer
	//color += vec4(offset,1)*0.2;

	return color;
}

struct sample
{
	vec4 color;
	float t;
};
sample accumulateSingleColorBlock_showOctree(sample s,BlockBox renderBlock,range r)
{
	//accumlate
	vec4 value, texValue;
	vec3 pos;	
	texValue = texture1D(tranf, renderBlock.singleColor);
	texValue.a=0.1;
	//analytical integrate
	float d=(r.tFar - s.t)/stepsize;
	float t = (1.0-texValue.a);
	float T = exp(-t*d);
	s.color.rgb += T*texValue.rgb;
	s.color.a += T;



	//safegard
	//int i=0;
	//while(s.color.a < 1.0 && s.t < r.tFar)
	//{
	//	if(i>200)
	//		break;
	//	i++;
	//	s.color.rgb += (1.0-s.color.a)*texValue.a*texValue.rgb;
	//	s.color.a += (1.0-s.color.a)*texValue.a;
	//	s.t += stepsize;
	//}

	return s;
}
 sample accumulateSingleColorBlock_R(sample s,BlockBox renderBlock,range r)
 {
	//accumlate
	vec4 value, texValue;
	vec3 pos;	
	texValue = texture1D(tranf, renderBlock.singleColor);

	//analytical integrate
	//float d=(r.tFar - s.t)/stepsize;
	//float t = (1.0-texValue.a);
	//float T = exp(-t*d);
	//s.color.rgb += T*texValue.rgb;
	//s.color.a += T;



	//safegard
	int i=0;
	while(s.color.a < 0.99 && s.t < r.tFar)
	{
		//if(i>500)
		//	break;
		i++;
		s.color.rgb += (1.0-s.color.a)*texValue.a*texValue.rgb;
		s.color.a += (1.0-s.color.a)*texValue.a;
		s.t += renderBlock.sampleRate;
	}

	return s;
 }

sample accumulateBlock_R1(vec3 rayO, vec3 rayDir, sample s,BlockBox renderBlock, range r)
{
	//warping the coordinate
	uvec3 blockPointer = renderBlock.blockPointer;

	float invboxsize =1.0/( renderBlock.boxMax.r - renderBlock.boxMin.r);
	vec3 offset = (1+1.0/(float(blockPoolDim*blockSize)-1))/float(blockPoolDim)*vec3(blockPointer);
	//vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	float shrinkFactor = (float(blockSize-1u)-1)/(float(blockSize-1u));
	//vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;

	//accumlate
	float tEnd = r.tFar;
	vec4 value, texValue;
	vec3 pos;	

	int i=0;
	while(s.color.a < 0.99 && s.t < tEnd)
	{
		//safeguard
		i++;
		if(i>200)
			break;

		pos = (rayO+s.t*rayDir);
		//return vec4(pos,1);
		
		pos = (pos-renderBlock.boxMin)*invboxsize;
		
		//return offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset);
		
		value = texture3D(blockPool, offset + (pos*invBlockPoolDim*shrinkFactor ));
		//value = texture3D(blockPool, offset + pos*shrinkFactor*invBlockPoolDim + shrinkOffset);
		texValue = texture1D(tranf, value.a);

		s.color.rgb += (1.0-s.color.a)*texValue.a*texValue.rgb;
		s.color.a += (1.0-s.color.a)*texValue.a;
		s.t += renderBlock.sampleRate;
		//s.color.a = 1;
		//return s;

	}

	//visualize block Pointer
	//color += vec4(offset,1)*0.2;

	return s;
}

sample accumulateBlock_R(vec3 rayO, vec3 rayDir, sample s,BlockBox renderBlock, range r)
{
	//warping the coordinate
	uvec3 blockPointer = renderBlock.blockPointer;

	float invboxsize =1.0/( renderBlock.boxMax.r - renderBlock.boxMin.r);
	vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	float shrinkFactor = (float(blockSize-1u))/(float(blockSize));
	vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;

	//accumlate
	float tEnd = r.tFar;
	vec4 value, texValue;
	vec3 pos;	

	int i=0;
	while(s.color.a < 0.99 && s.t < tEnd)
	{
		//safeguard
		i++;
		if(i>200)
			break;

		pos = (rayO+s.t*rayDir);
		//return vec4(pos,1);
		
		pos = (pos-renderBlock.boxMin)*invboxsize;
		
		//return offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset);
		
		//value = texture3D(blockPool, offset + (pos*invBlockPoolDim*shrinkFactor ));
		value = texture3D(blockPool, offset + pos*shrinkFactor*invBlockPoolDim + shrinkOffset);
		texValue = texture1D(tranf, value.a);

		s.color.rgb += (1.0-s.color.a)*texValue.a*texValue.rgb;
		s.color.a += (1.0-s.color.a)*texValue.a;
		s.t += renderBlock.sampleRate;
		//s.color.a = 1;
		//return s;

	}

	//visualize block Pointer
	//color += vec4(offset,1)*0.2;

	return s;
}
vec4  octreeRestartRayCasting(vec3 start, vec3 backPosition)
{
	vec4 color=vec4(0.0,0.0,0.0,0.0);
	vec3 rayDir = normalize(backPosition - start);
	float tMax = length(backPosition - start);
	float t = 0;//15*0.0014;
	
	int i=0;
	vec3 rayO = start;
	BlockBox renderBlock;
	range r ;
	vec3 pos;
	//sample s;
	//s.t = t;
	//s.color = vec4(0,0,0,0);
	while(t<tMax-0.01 && color.a <0.99)
	{
		//ray decending, find the bounding box for current block
		pos = rayO + t*rayDir;
		renderBlock = findCurrentBlockAndBoundingBox(pos+0.0001*rayDir);
		//calculate ray box intersection

		r = intersectBox(rayO , rayDir, renderBlock.boxMin, renderBlock.boxMax);

		if(i==30)
		{
			//return vec4(1,0,0,1);
			break;
		}
		i++;

		//accumulate volume
		if(renderBlock.hasBlock)
		{
			//color = vec4(rayO+r.tNear*rayDir,1);
			color = accumulateBlock(rayO, rayDir, color,renderBlock, r);
			//s = accumulateBlock_R(rayO, rayDir, s,renderBlock, r);
		}
		else
		{
			color = accumulateSingleColorBlock(color,renderBlock,r);
		}

		//advance ray
		//t = s.t;
		t = r.tFar;
	}

	return color;
}

vec4  octreeRestartRayCasting_R(vec3 start, vec3 backPosition)
{
	vec3 rayDir = normalize(backPosition - start);
	float tMax = length(backPosition - start);
	
	int i=0;
	vec3 rayO = start;
	BlockBox renderBlock;
	range r ;
	vec3 pos;

	sample s;
	s.t = 15*0.0014;//stepsize;
	s.color = vec4(0,0,0,0);

	while(s.t<tMax-0.01 && s.color.a <0.99)
	{
		//ray decending, find the bounding box for current block
		pos = rayO + s.t*rayDir;
		renderBlock = findCurrentBlockAndBoundingBox(pos);
		//calculate ray box intersection

		r = intersectBox(rayO , rayDir, renderBlock.boxMin, renderBlock.boxMax);

		if(i==40)
		{
			//return vec4(1,0,0,1);
			break;
		}
		i++;

		//accumulate volume
		if(renderBlock.hasBlock)
		{
			//color = vec4(rayO+r.tNear*rayDir,1);
			//color = accumulateBlock(rayO, rayDir, color,renderBlock, r);
			s = accumulateBlock_R(rayO, rayDir, s,renderBlock, r);
		}
		else
		{
			s = accumulateSingleColorBlock_R(s,renderBlock,r);
		}

	}

	return s.color;
}
/////////////////////////////////////////////////////////////////////////////////////

vec4 octreeRestartRayCasting_R1(ray r, range T)//r.o ,r.d
{
	vec4 color = vec4(0,0,0,0);
	float t = T.tNear;
	float tMax = T.tFar;
	vec3 pos;
	range tRange;
	BlockBox renderBlock;

	while(t<tMax-0.001 && color.a <0.99)
	{
		pos = r.o + t*r.d;
		pos = (pos+vec3(1,1,1))*0.5;
		return vec4(pos,1);

		renderBlock = findCurrentBlockAndBoundingBox(pos+0.00001*r.d);
		//calculate ray box intersection

		tRange = intersectBox(r.o , r.d, renderBlock.boxMin, renderBlock.boxMax);


	}

	return color;	
}






///////////////////////////////////////////////////////////////////////////////////////
//===================================================================
vec3 OctreeCoord(vec3 pos)
{ 
	uvec3 nodeClusterPointer=uvec3(0u,0u,0u);
	uvec3 childrenIndex; //@@@@you can't write childrenIndex = uvec3(2*pos) !!!!!!!!!!WTF@@@@

	nodeTexel currentNode;
	currentNode.first = 1u<<31u; currentNode.second = 0u;

	uint octreeLevel = 0u;
	//float scaleFactor = 1.0; //@@ the concept of scale is wrong @@!!!!!!
	
	bool SubD=true;
	int c=0;
	while(SubD) //@@can't write isSubDividable(currentNode);@@!!!WTF
	{
		c++;
		if(c==10)
			break;
		//decending
		childrenIndex.x = uint((2*pos.x));
		childrenIndex.y = uint((2*pos.y));
		childrenIndex.z = uint((2*pos.z));
		pos = fract(2*pos);
		
		currentNode=getChildrenNode( getNodeClusterPointer(currentNode), childrenIndex.z, childrenIndex.y, childrenIndex.x);
		SubD = isSubDividable( currentNode );

		//scaleFactor *= 0.5;
	}
	//calculate blockoffset in the pool
	uvec3 blockPointer=getBlockPointer(currentNode);

	if(!hasBlock(currentNode))
		return vec3(0.0,0.0,0.0);

	//boundary shrink factor
	vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	float shrinkFactor = (float(blockSize)-1)/(float(blockSize));
	vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;
	//scale the pos
	return offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset);
	//return offset + pos*invBlockPoolDim;
}

vec4 DirectVolumeRenderingOctree(vec3 front, vec3 back, float upperThreshold, float threshold, float offset)
{
	vec3 dir;
	dir = back - front;
	float len = length(dir);
	vec3 norm_dir = normalize(dir);
	//float iso_value = 0.1;

	vec3 rayIncrement = stepsize*norm_dir;

	vec3 pos = front;
	vec4 value=vec4(0,0,0,0);
	// accumlate value
	vec4 result = vec4(0,0,0,0);//RGBA
	float raylength = 0;
//####
	//return vec4(OctreeCoord(pos+norm_dir*0.1),1.0);
//####

    float scale = 1.0/(upperThreshold - threshold);
    float res = 0;
    
	for(int i=0; i<3000; i++)
	{
		pos += rayIncrement;

		//value = texture3D(blockPool, pos);
		//if(pos.y>isoValue*1.6)
			value = texture3D(blockPool, OctreeCoord(pos));
		//else
		//	value = vec4(0,0,0,0);
	
		raylength += length(rayIncrement);

			vec4 texValue = texture1D(tranf, value.a);

			result.rgb += (1.0-result.a)*texValue.a*texValue.rgb;
			result.a += (1.0-result.a)*texValue.a;

			if (result.a >= 0.99 || raylength > len-0.01 ) 
			{
				break;
			}
	}
	return result;

}
//===================================================================
vec4 visualizeOctree(vec3 start, vec3 end, vec3 O, vec3 screenNorm)
{
	vec4 color;
	float viewt = (1.0/dot(screenNorm,normalize(start-O)));
	float t = length(start-O);
	//viewt/t

	color = vec4(vec3(viewt/t),1);
	//return vec4(screenNorm,1);
///////////////////////////////////////////////////////
	vec3 rayDir = normalize(end - start);
	float tMax = length(end - start);
	
	int i=0;
	vec3 rayO = start;
	BlockBox renderBlock;
	range r ;
	vec3 pos;

	sample s;
	s.t = 15*0.0014;//stepsize;
	s.color = vec4(0,0,0,0);

	while(s.t<tMax-0.01 && s.color.a <0.99)
	{
		//ray decending, find the bounding box for current block
		pos = rayO + s.t*rayDir;
		renderBlock = findCurrentBlockAndBoundingBox_MultiRes(pos, viewt, t+s.t, O);
		//calculate ray box intersection

		r = intersectBox(rayO , rayDir, renderBlock.boxMin, renderBlock.boxMax);

		if(i==40)
		{
			//return vec4(1,0,0,1);
			break;
		}
		i++;

		//accumulate volume
		if(renderBlock.hasBlock)
		{
			//color = vec4(rayO+r.tNear*rayDir,1);
			//color = accumulateBlock(rayO, rayDir, color,renderBlock, r);
			//s = accumulateBlock_R(rayO, rayDir, s,renderBlock, r);
			s = accumulateSingleColorBlock_showOctree(s,renderBlock,r);
		}
		else
		{
			s = accumulateSingleColorBlock_showOctree(s,renderBlock,r);
		}

	}

	return s.color;
	//return color;
}
//////////////////////////////////////////////////////////////////

vec4 multiResOctree(vec3 start, vec3 end, vec3 O, vec3 screenNorm)
{
	vec4 color;
	float viewt = (1.0/dot(screenNorm,normalize(start-O)));
	float t = length(start-O);
	//viewt/t

	color = vec4(vec3(viewt/t),1);
	//return vec4(screenNorm,1);
///////////////////////////////////////////////////////
	vec3 rayDir = normalize(end - start);
	float tMax = length(end - start);
	
	int i=0;
	vec3 rayO = start;
	BlockBox renderBlock;
	range r ;
	vec3 pos;

	sample s;
	s.t = 15*0.0014;//stepsize;
	s.color = vec4(0,0,0,0);

	while(s.t<tMax-0.01 && s.color.a <0.99)
	{
		//ray decending, find the bounding box for current block
		pos = rayO + s.t*rayDir;
		renderBlock = findCurrentBlockAndBoundingBox_MultiRes(pos, viewt, t+s.t, O);
		//calculate ray box intersection

		r = intersectBox(rayO , rayDir, renderBlock.boxMin, renderBlock.boxMax);

		if(i==40)
		{
			//return vec4(1,0,0,1);
			break;
		}
		i++;

		//accumulate volume
		if(renderBlock.hasBlock)
		{
			//color = vec4(rayO+r.tNear*rayDir,1);
			//color = accumulateBlock(rayO, rayDir, color,renderBlock, r);
			s = accumulateBlock_R(rayO, rayDir, s,renderBlock, r);
			//s = accumulateSingleColorBlock_showOctree(s,renderBlock,r);
		}
		else
		{
			s = accumulateSingleColorBlock_R(s,renderBlock,r);
			//s = accumulateSingleColorBlock_showOctree(s,renderBlock,r);
		}

	}

	return s.color;
	//return color;
}

void main()
{
	vec3 boxMin = vec3(-1.0, -1.0, -1.0);
    //vec3 boxMin = vec3(0, 0, 0);
    vec3 boxMax = vec3(1.0, 1.0, 1.0);

	vec4 O = (invViewMatrix*vec4(0,0,0,1));
	vec3 eyeO = O.xyz/O.w;
	mat3 rotateMat = mat3(vec3(invViewMatrix[0]),vec3(invViewMatrix[1]),vec3(invViewMatrix[2]));
	vec3 eyeD = rotateMat*normalize(vec3(uv,-2.0));

	vec3 screenNorm = rotateMat*normalize(vec3(0,0,-2.0));
	
	range t = intersectBox(eyeO, eyeD, boxMin, boxMax);

	vec3 start = vec3(0);
	vec3 end = vec3(0);
	
	if(t.tNear<0)
		t.tNear = 0;
	
	if(t.tNear>t.tFar)
		gl_FragColor = vec4(0,0,0,1);
	else
	{
		//gl_FragColor = vec4(pos*0.5+0.5,t.tNear);
		start = (eyeO + eyeD*t.tNear)*0.5+0.5;
		end = (eyeO + eyeD*t.tFar)*0.5+0.5;
		vec3 correctedO = eyeO*0.5+0.5;
		//gl_FragColor = vec4(end,1);


		//gl_FragColor = vec4(rotateMat[2],1);	
		//gl_FragColor = vec4(start,1);
		//gl_FragColor = DirectVolumeRendering(start, end, 0, 0, 0);
		ray r;
		r.o = eyeO; r.d = eyeD;
		//gl_FragColor = octreeRestartRayCasting_R1(r, t);//r.o ,r.d
		if(isoValue > 0.5)
			gl_FragColor = visualizeOctree(start, end, correctedO,screenNorm);
			//gl_FragColor = octreeRestartRayCasting_R(start, end);
			//gl_FragColor = octreeRestartRayCasting(start, end);
		else
			gl_FragColor = multiResOctree(start, end, correctedO,screenNorm);
		//gl_FragColor = invViewMatrix[3];
	}
}